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We present a study of hydrogen at pressures higher than molecular dissociation using the Coupled 
Electron-Ion Monte Carlo method. These calculations use the accurate Reptation Quantum Monte 
Carlo method to estimate the electronic energy and pressure while doing a Monte Carlo simulation 
of the protons. In addition to presenting simulation results for the equation of state over a large 
region of phase space, we report the free energy obtained by thermodynamic integration. We find 
very good agreement with DFT calculations for pressures beyond 600 GPa and densities above 
p = lAg/cm 3 . Both thermodynamic as well as structural properties are accurately reproduced by 
DFT calculations. This agreement gives a strong support to the different approximations employed 
in DFT, specifically the approximate exchange-correlation potential and the use of pseudopotentials 
for the range of densities considered. We find disagreement with chemical models, which suggests a 
reinvestigation of planetary models, previously constructed using the Saumon-Chabrier-Van Horn 
equations of state. 



I. INTRODUCTION 

Although hydrogen is the first clement in the periodic 
table, its phase diagram has diverse phases. As the most 
abundant element in the universe, it is important to have 
an accurate understanding of its properties for a large 
range of pressure and temperature. A qualitative de- 
scription is not sufficient because; for example, models 
of hydrogenic planets require accurate results to make 
correct predictions pQ. 

The high pressure phases of hydrogen have received 
considerable attention in recent years, both from the- 
ory and experiment. At lower temperatures, static com- 
pression experiments using diamond anvil cells can reach 
pressures of 320 GPa, where the quest to find the metal- 
insulator and molecular-atomic transitions in the solid 
phase still continues |2J. Dynamic compression experi- 
ments using either isentropic compression or shock waves, 
are used at higher temperatures and can now reach pres- 
sures above 200 GPa [31 0] . Even though experimental 
techniques at high pressure have improved considerably 
over the last decade, they are still not accurate enough 
to provide conclusive answers to many of the relevant 
questions. Although this situation might change in the 
near future with the construction of more powerful ma- 
chines such as the National Ignition Facility, computer 
simulations today provide the most reliable method for 
determining the thermodynamic properties at high pres- 
sures and temperatures. 
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Many theoretical techniques have been used includ- 
ing: free energy minimization methods in the chemical 
picture [5J E], restricted Path Integral Monte Carlo 
(PIMC) [8] and density functional theory (DFT) based 
molecular dynamics (MD) [gj [TOl ED EH ECU H EES] . All 
of these methods employ different approximations that 
can affect properties in ways that are difficult to quan- 
tify due to the lack of conclusive experimental results. 
While free energy methods are typically accurate in the 
molecular phase at low pressures, where molecules are 
tightly bound and there are enough experimental results 
to produce accurate empirical potentials, at higher den- 
sity, with the onset of dissociation and metallization in 
the liquid, they become unreliable. Restricted PIMC, on 
the other hand, is accurate at very high temperatures 
where the nodes of the density matrix are known, but 
at temperatures below approximately 20,000K its accu- 
racy (and efficiency) has been limited. For intermedi- 
ate temperatures and high pressures, DFT has become 
the computational method of choice over the last decade, 
mainly due to its advanced development stage and easy 
accessibility with many available codes. Practical im- 
plementations of DFT employ pseudopotentials and ap- 
proximate exchange-correlation functionals and do not 
typically estimate quantum proton effects; these approx- 
imations limit its accuracy and applicability, especially 
at high pressures. Despite its possible limitations, DFT 
is state of the art in ab initio simulations. For exam- 
ple, planetary models are being built with its equation of 
state, superseding the well known Saumon-Chabrier-Van 
Horn (SCVH) multiphase equation of state [5l [16]. 

Because of its wide spread use and potential impact 
in the near future, it is important to test the validity of 
the DFT approximations at extreme conditions, and to 
determine its range of applicability. In order to obtain 
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more accurate results, especially at intermediate temper- 
atures, we need to employ methods that can go beyond 
the usual single-body mean field approximations typi- 
cally used. Quantum Monte Carlo (QMC) is a perfect 
candidate for the task, presenting a good balanced be- 
tween speed and accuracy. It does not rely on pseudopo- 
tentials, and correlation effects are treated explicitly in 
the full many-body problem. At present, QMC has the 
potential to be more accurate for electronic properties 
than DFT while being considerably less expensive than 
quantum chemistry methods [T7J HH1 [TH] . 

Coupled Electron-Ion Monte Carlo (CEIMC) is a QMC 
based ab initio method developed to use QMC electronic 
energy in a Monte Carlo simulation of the ionic degrees 
of freedom [201 121] ■ Thanks to recent advances in QMC 
methodology we can now obtain results with small sys- 
tematic errors. Specifically, the use of Twist Averaged 
Boundary Conditions (TABC) [23] together with recently 
developed finite-size correction schemes [53] |25] allow us 
to produce energies that are well converged to the ther- 
modynamic limit with ~100 atoms; see appendix |B] for 
additional details. Improvements in the wavefunctions 
used for high pressure hydrogen allow us to get very ac- 
curate results |26j and avoid the main limitations of pre- 
vious applications of the CEIMC method. 

The impressive increase in computational power in re- 
cent years made possible a comprehensive study of the 
equation of state of hydrogen. In this paper we present 
the results of free energy calculations of hydrogen in its 
atomic liquid phase for pressures beyond 150 GPa. In 
section [TTJ we present a brief description of the CEIMC 
method. In section |III| and appendix [X] we present the 
equation of state (EOS), including the free energy cal- 
culations. In section |IV| we present a comparison with 
other methods, with a special attention to DFT-based 



V] we draw some conclu- 
andOwe describe details 



MD results, while in section 
sions. Finally in appendices [B 
of the method to estimate the finite size corrections and 
of the RQMC implementation, respectively. 



II. COUPLED ELECTRON-ION MONTE 
CARLO METHOD 



CEIMC, in common with the large majority of ab initio 
methods, is based on the Born-Oppenheimer separation 
of electronic and ionic degrees of freedom. In addition, 
the electrons are considered to be in their ground state, 
for some particular arrangement of the protons. Protons, 
either considered as classical or quantum particles [37], 
are instead assumed to be at thermal equilibrium with a 
heat bath at a temperature T. In the present calculation, 
the system of N protons and N electrons is enclosed in 
a fixed volume V at a number density n = N/V, which 
we often express with the parameter r s = (3/47rn)^ 1/ ' 3 -'. 
The mass density is related to r s by: p = (3m^)/(47rr^), 
with mh the mass of a hydrogen atom. 



We start from the non-relativistic Hamiltonian 

2N „2 
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Zi, mi, fi represent respectively the valence, mass and 
position operators of particle i, and Aj — ti 2 /2m i . Let us 
denote with R — [fx, ■ ■ ■ , rjv) and S — (rjv+i, • • ■ , ?2n) 
the set of coordinates of all electrons and protons respec- 
tively. [37] 

Within the Born-Oppenheimer approximation, the en- 
ergy of the system for a given nuclear state S is the expec- 
tation value of the hamiltonian H over the corresponding 
exact ground state | < &o(S')) 



(2) 



Ebo(S) = ($ (S) H *o(5)>, 



which is a 3./V-dimensional integral over the electron co- 
ordinates in configuration space 



Ebo(S) 



dR $* (R\S)H(R, S)<S> (R\S) 
dR\<f> (R\S)\ 2 E L (R\S), 



with the local energy defined as 



E L (R\S) = 



H(R, S)<S> (R\S) 
$a(R\S) 



(3) 



(4) 



In this work, we use Reptation Quantum Monte Carlo 
(RQMC) [28] with the bounce algorithm [29] to solve the 
electronic problem. 

With the ability to compute the Born-Oppenheimer 
electronic energy, the Metropolis algorithm is able to gen- 
erate a sequence of ionic states according to the Boltz- 
mann distribution P(S) cx e~P E Bo(S) a ^ ^ ne mverse tem- 
perature p. In CEIMC the estimate of Ebo(S) for a 
given trial function is computed by QMC and it is there- 
fore affected by statistical noise which, if ignored, will 
bias the result. In the Penalty Method [22 we require de- 
tailed balance to hold on average (over the noise distri- 
bution). The noise on the energy difference causes extra 
rejection with respect to the noiseless case. 

The accuracy of the calculations depend crucially on 
the choice of wavefunction. The Slater- J astrow wave- 
function has the form: 



(5) 



where U is the sum over all distinct pairs of particles 
of the RPA-Jastrow function [30] • The orbitals in the 
Slater determinant are obtained from a DFT calculation 
in a planewave basis, using the local density approxima- 
tion for the exchange-correlation functional, as parame- 
terized by Perdew, et. al. [31] [32]. The resulting orbitals 
are transformed to a cubic spline basis from which they 
are interpolated in the QMC calculations. The use of 
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FIG. 1: Free energy surface of hydrogen as a function of 
temperature and density. 



a spline basis in the RQMC runs represents a large in- 
crease in the efficiency of the simulations. To reduce the 
computational overhead produced by the DFT calcula- 
tions, which must be converged to self-consistency for 
every ionic step, we perform a single self-consistent DFT 
calculation at the gamma point. Using the resulting elec- 
tron density, we build and diagonalize the Kohn-Sham 
hamiltonian at the specific points in the Brillouin zone 
used in the TABC calculations. 

We apply a backflow transformation to the electron co- 
ordinates in the Slater determinant; this introduces cor- 
relations and improves the nodal surfaces of the DFT 
orbitals. The form of the transformation is: 

Si = fi + ^2 rjij ( i m \)fij, (6) 

3 

where rj is cither a parameterized electron-electron back- 
flow function or an analytic forms derived using Bohm- 
Pines collective coordinates approach [33] . As shown in a 
recent publications [T5J|2B], the use of backflow transfor- 
mations improves the trial wavefunction, especially for 
small projection times. The use of DFT orbitals with 
RPA-derived Jastrow and backflow functions represent a 
good balance between accuracy and efficiency. 



r s AE (mHa) AP (GPa) 
TfJ5 4.0 (7) 7 (3) 
1.10 3.8 (3) 9 (1) 
1.25 1 2.8 (5) | 5 (1) 

TABLE I: Corrections to the energy and pressure of hydro- 
gen from quantum effects of the protons, from PIMC sim- 
ulations with CEIMC, at a temperature of 2000 K: AE = 
(E - E c iassicai)/N and AP = P — Pciassical- Errors in paren- 
theses. 



length used in RQMC were chosen to reach a convergence 
of the energy of 0.2-0.3 mHa/atom; see appendix [C| for 
details. We used TABC with a grid of 64 twists (96 for 
r s < 1.10), which together with the use of recently de- 
veloped finite-size correction schemes for QMC energies, 
allows a significant reduction of size effects; see appendix 
|B]for details and discussion. 

We performed a 36 CEIMC simulations (not including 
those related to the coupling constant integration), the 
results are reported in Table III of appendix [AJ Our sim- 
ulations were first equilibrated using effective pair poten- 
tials built from reflected Yukawa functions, which were 
chosen such they reproduced the radial distribution func- 
tions of the QMC systems. We then performed 2000-3000 
equilibration steps with CEIMC. Statistics were gathered 
in the following 5000-15000 steps, the number depending 
on temperature and density. 

The protons are treated as classical particles in the re- 
sults presented in Table III and in the free energy calcu- 
lations discussed below. Quantum effects of the protons 
could be important at low temperatures at the densities 
considered in this work. In order to assess their effect on 
the thermodynamic properties, we performed Path Inte- 
gral Monte Carlo (PIMC) calculations for the protons on 
the potential energy surface defined by the zero temper- 
ature RQMC method. This is an extension of CEIMC to 
path integral calculations [2T]. The resulting corrections 
to the energy and pressure are given in Table I, at T= 
2000K for 3 densities. At this temperature, which is the 
lowest temperature studied in this work, the corrections 
to the pressure are below 1%. 



III. EQUATION OF STATE OF HYDROGEN 

In this paper we calculate the free energy of hydrogen 
as a function of temperature and density. We explore 
the range 2 000K < T < 10 000K and 0.7 g cm" 3 < p < 
2.4 g cnr 3 . For p < 0.7 g cm~ 3 , molecular dissociation 
becomes the dominant feature in the EOS, requiring a 
more detailed study than the one performed here to reach 
a similar level of accuracy. Simulations of the molecular 
liquid are in progress. 

The CEIMC calculations were performed with 54 hy- 
drogen atoms using RQMC. The time step and projection 



A. Free Energy Integration 

We used Coupling Constant Integration (CCI) [M] to 
calculate the free energy of hydrogen at a reference point 
chosen as T=6000K and r s — 1.25. In CCI a A-dependent 
potential energy, U(A), is a linear combination of two dif- 
ferent potentials: V(X) — Vo + A(Vi — Vq). The difference 
in free energy between systems and 1 is then: 

Fi-F = f dX((yx-V Q )) vw . (7) 
Jo 
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FIG. 2: The entropy per atom as a function of temperature 
and density. 



We chose the reference potential to be a reflected Yukawa 
pair potential: 



Vo(r) 



-b(L-r) -bL/2 

r<L/2, 







r > L/2, 



(8) 



where b=2.5 a.u. and L is the length of the simula- 
tion cell. The reflection makes the function and its first 
derivative continuous at the cut-off, r = L/2. We first 
used CCI to calculate the free energy difference between 
the Yukawa potential and a system of non-interacting 
particles. We performed a second CCI with CEIMC to 
calculate the free energy difference between the Yukawa 
model and the QMC system. We obtained a value of - 
0.5737(1) Ha/atom for the free energy and 1.98(1)*1(T 5 
(Ha/K) for the entropy at the reference point. 



is much smaller than this. Figures [T] and [2] show plots 
of the free energy and entropy in the region of the phase 
diagram studied. The coefficients of the expansion are 
given in Table II. 
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FIG. 3: Comparison of EOS data with DFT-based BOMD 
simulations. 



IV. COMPARISON WITH OTHER METHODS 



B. Fits to EOS 

The free energy as a function of temperature and den- 
sity was obtained by using the functional form: 



F(T,p) = ^2^2 c tk 9i(T) g k (p), 



(9) 



i=i k=i 



where gi{x) — {1, x, x 2 , xlnx, x 2 Inx}. We do a least 
squares fit between the analytical function and the 
CEIMC data using derivatives of this expansion: 



P(T,P) 
E(T,p) 
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(10) 

(11) 



The free energy fit reproduces the energies and pres- 
sures obtained from the CEIMC simulations to within 
0.5% and 1.5% respectively, although the average error 



Figure [3] shows a comparison of the pressure and en- 
ergy between CEIMC simulations and DFT based Born- 
Oppenheimer Molecular Dynamics (BOMD) simulations. 
The BOMD simulations were performed in the NVT- 
ensemble (with a weakly coupled Berendsen thermostat) 
using the Qbox code [35]. We used the Perdew-Burke- 
Ernzerhof (PBE) exchange-correlation functional and a 
Hamman type [35] local pseudopotential with a core ra- 
dius of r c = 0.3 a.u. to represent hydrogen. The simu- 
lations were performed with 250 hydrogen atoms using 
a plane-wave cutoff of 90 Ry (115 Ry for rs > 1.10)with 
periodic boundary conditions (r point). Corrections to 
the EOS were added to extrapolate results to infinite cut- 
off and to account for the Brillouin zone integration. To 
do this we studied 15-20 statistically independent static 
configurations at each density by using a 4x4x4 grid of 
k-points with a plane-wave cutoff of at least 300 Ry. See 
ref. [36] for additional details of the DFT simulations. 

There is a good agreement between the two meth- 
ods, especially at higher densities where the difference 
in pressure is within error bars. At lower densities, the 
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TABLE II: Coefficients of the expansion of the free energy; energy in Hartree/atom, temperature in K and density in g cm 



pressure difference increases reaching an average of ap- 
proximately 5% close to the dissociation regime. There 
is less reason to expect exact agreement of the energies 
calculated since the DFT calculations use pseudopoten- 
tials as well as approximate exchange-correlation func- 
tional which can modify the zero of energy. However, 
the temperature and density dependence is well repro- 
duced with an almost uniform energy difference of 0.8% 
on the entire phase diagram. Figure [4] shows a compar- 
ison of the proton-proton distribution function for sev- 
eral thermodynamic conditions. The agreement between 
the two methods is again remarkable. The structure of 
the liquid is reproduced by the DFT simulations quite 
accurately, even the short range correlation peak that 
develops at the lower temperatures and higher densities. 
Figure [5] shows a comparison of the entropy as a function 
of density along several isotherms. We can see that for 
densities beyond p — lAg/cm 3 , the entropies obtained 
with the two methods are undistinguishablc. In general, 
we obtain very good agreement between the two methods 
for pressures beyond 600 GPa. At lower pressures, the 
agreement is not perfect, but still very good, with BOMD 
predicting a slightly higher entropy than CEIMC. We are 
currently expanding our calculations to lower density to 
provide an additional benchmark of the DFT method 
close to the molecular dissociation region. 

Figure [6] shows a comparison of the pressure as a func- 
tion of density obtained with CEIMC and BOMC and the 
SCVH equation of state at T=6000 K. At pressures well 
below the dissociation regime, SCVH produces very good 
results, but at higher densities the model can not capture 
all the features which result from strong atomic coupling; 
it predicts a qualitatively different behavior with higher 
pressures at lower densities. 
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FIG. 4: Comparison of radial distribution functions between 
BOMD (blue) and CEIMC (red). 



V. CONCLUSIONS 

In summary, we performed a comprehensive study of 
the free energy and the equation of state of warm dense 
liquid hydrogen in its atomic phase using energies com- 
puted with quantum Monte Carlo methods. We provide 
a fit to the free energy which can be used as input to 
models of Jovian planets or in the formulation of more 
accurate chemical models. Given the current status of 
DFT and its possible limitations at these conditions, it 
is crucial to benchmark its predictions against more ac- 



curate methods. We provide such a critical test. Our 
results indicate that DFT-based BOMD simulations pro- 
vide a very good description of both thermodynamic and 
structural properties of hydrogen for the studied condi- 
tions. The equation of state of SCVH, used in the study 
of planetary interiors for more than a decade, is shown 
to produce inaccurate results in the atomic regime. This 
suggests that planetary models should be reinvestigated 
with a more accurate equations of state, such as the one 
presented in this work. 
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tion of density computed with SCVH, BOMD and CEIMC at 
T=6000 K. 



APPENDIX A: EQUATION OF STATE TABLE 

Table III gives the energy per atom and pressure as cal- 
culated by CEIMC on a grid in temperature and density. 
Finite size corrections (sec next section) arc included in 
the reported results, as well as corrections to the pres- 
sure from non-zero time step and finite projection time in 
RQMC. The latter corrections are negligible in the case 
of total energies. 



T(K) 


r s 


Energy (Ha) P 


ressure 


(GPa) 


2000 


1.05 


-0.3846 


(3) 


1576 


(2) 


3000 


1.05 


-0.3777 


(2) 


1607 


(1) 


4000 


1.05 


-0.3707 


(5) 


1640 


(3) 


6000 


1.05 


-0.3569 


(4) 


1701 


(2) 


8000 


1.05 


-0.3458 


(4) 


1753 


(2) 


10000 


1.05 


-0.3316 


(8) 


1814 


(4) 


2000 


1.10 


-0.4170 


(2) 


1157 


(2) 


3000 


1.10 


-0.4097 


(2) 


1190 


(1) 


4000 


1.10 


-0.4026 


(3) 


1219 


(1) 


6000 


1.10 


-0.3898 


(4) 


1270 


(2) 


8000 


1.10 


-0.3777 


(5) 


1315 


(2) 


10000 


1.10 


-0.3660 


(4) 


1362 


(2) 


2000 


1.15 


-0.4419 


(2) 


861 


:i) 


3000 


1.15 


-0.4356 


(3) 


883 


;2) 


4000 


1.15 


-0.4285 


(7) 


911 


;s) 


6000 


1.15 


-0.4151 


(6) 


956 


;3) 


8000 


1.15 


-0.4018 


(9) 


1003 


(3) 


10000 


1.15 


-0.3891 


(8) 


1048 


(2) 


2000 


1.25 


-0.4790 


(2) 


485 


:i) 


3000 


1.25 


-0.4721 


(2) 


504 


;i) 


4000 


1.25 


-0.4673 


(4) 


517 


2) 


6000 


1.25 


-0.4549 


(6) 


556 


;i) 


8000 


1.25 


-0.4419 


(7) 


590 


:s) 


10000 


1.25 


-0.4324 


(7) 


619 


:2) 


2000 


1.40 


-0.5117 


(4) 


214 


;i) 


3000 


1.40 


-0.5057 


(2) 


222 


;i) 


4000 


1.40 


-0.4993 


(5) 


234 


;i) 


6000 


1.40 


-0.4869 


(5) 


257 


3) 


8000 


1.40 


-0.4767 


(5) 


277 


;i) 


10000 


1.40 


-0.4674 


(4) 


297 


:i) 


2000 


1.55 


-0.5330 


(2) 


111 


:i) 


3000 


1.55 


-0.5230 


(2) 


105 


;i) 


4000 


1.55 


-0.5157 


(3) 


117 


;i) 


6000 


1.55 


-0.5027 


(4) 


134 


;i) 


8000 


1.55 


-0.4938 


(2) 


143 


;i) 


10000 


1.55 


-0.4831 


(6) 


163 


;2) 



TABLE III: Energy/atom and pressure as calculated with 
CEIMC. Statistical errors are given in parentheses. 



APPENDIX B: FINITE SIZE EFFECTS 

Due to the high computational demands of QMC, our 
simulations are restricted to systems with 128 atoms or 
fewer. Many techniques have been developed in order to 
obtain useful results with finite systems. In this work we 
use TABC (the generalization of Brillouin zone integra- 
tion to many-body quantum systems) to eliminate shell 
effects in the kinetic energy of metallic systems. Twisted 
boundary conditions when an electron wraps around the 



simulation box are defined by: 

*<,[■->■] ■ a....! -"'>!',(..../•;....}. (bi) 

Observables are then averaged over the all twist vectors, 
similar to one-body theories: 

<A>= £ j^<* e \A\* e >. (B2) 



7 



This has been shown to restore the 1/N dependence of 
the energy in QMC calculations, absent when PBC are 
used. 

As first shown by Chiesa et.al. [M], most of the re- 
maining finite size errors in the potential and kinetic en- 
ergies of QMC simulations come from discretization er- 
rors induced by the use of PBC. To see this, notice that 
we can write the potential energy of a system of N elec- 
trons as: 



(B3) 



where f2 is the volume of the supercell, v(k) is the 
Fourier transform of the Coulomb potential, S^^k) = 
< p(k)p(~k) > is the structure factor, p(k) — JV Zie lk ' ri 

[3"5] and k is the set of lattice vectors in reciprocal 
space of the supercell. As we approach the thermo- 
dynamic limit N — > oo, the structure factor converges 
(Sj^(k) — > Soo(k)) and the sum becomes an integral: 

V Sfc^o ~~ * I (27r) 3 • If we assume that the structure 
factor is essentially converged for some finite number of 
atoms, then most of the finite size error in the potential 
energy comes from the omission of the k = term in the 
sum. This can be estimated using the Poisson summation 
formula: 



(B4) 



We know from the RPA, exact in the limit of k — > 0, that 
S(k) ss k 2 as k — > 0. The leading order correction to the 
potential energy is: 
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By estimating the structure factor for small k during our 
simulation and extrapolating its behavior to k=0, we ob- 
tain the desired corrections. In the case of the kinetic 
energy, following a similar argument we obtain for the 
correction to second order 1^1 1251: 
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where u(k) is the electron-electron Jastrow; we used its 
RPA form. 

To check the finite-size corrections for dense hydrogen, 
we performed simulations with 32, 54 and 108 atoms at 
r s =1.85 and r s =1.25. We used TABC with Variational 
Monte Carlo energies with 108 twists for the 32 atom 



system and 32 twists for the 54 and 108 atoms. Figure [7] 
shows a comparison of the radial distribution functions 
between the systems with different number of atoms, for 
the two densities studied. As can be seen, the agree- 
ment between the 3 simulations is very good, with no 
noticeable difference between the systems with 54 and 
108 atoms. Using the fact that TABC restores the 1/N 
dependence of the properties, we can compare the results 
of the size correction formulas with a 1/N extrapolation. 
Table IV shows a comparison of finite size corrections 
taking the system with 54 atoms as the reference. As 
can be seen, the correction for the lower density system 
agrees very well with the extrapolated value. In the case 
of the higher density system, the agreement is less good 
but still acceptable. We attribute the disagreement at 
higher densities to the differences in TABC used for the 
systems with different sizes. 



r s 


Energy (mHa) 


Pressure (GPa) 




AE N 


AE S 


AP N 


AP S 


1.25 


7.7 


10.0 


10.8 


17.6 


1.85 


5.4 


5.5 


3.3 


3.1 



TABLE IV: Comparison of finite-size corrections between 
a size extrapolation (AEn ,APn) and formulas B5 and B6 
(AEsAPs)- 



APPENDIX C: DETAILS OF RQMC 
CALCULATIONS 



As mentioned in the text, the parameters of the 
RQMC calculations were varied with density to obtain 
uniform accuracy over the whole phase diagram. Table 
V shows the time steps and projection times used in the 
simulations for the different densities studied. As is well 
known, total energies converge faster as a function of 
imaginary time because their convergence error is second 
order with respect to the trial wavefunction. The kinetic 
and potential energies, and hence the pressure, are only 
first order, which means that longer projection times are 
needed to obtain converged results. During the course 
of the simulations, we only require accurate energies, so 
a smaller projection time is used; one that is insufficient 
for accurate pressures. In order to obtain converged 
results for the pressure, we calculated a correction for 
the finite projection time and non-zero time step by 
studying approximately a dozen protonic configurations 
at each density. The corrections were found to be 
independent of the precise proton configuration but 
density dependent. The corrections to the total energy 
are negligible within errors. 
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FIG. 7: Comparison of proton-proton radial distribution 
functions between systems with different number of electrons. 
In the upper panel, the temperature is 4000K and r s =1.25. In 
the lower panel, the temperature is 3000K and r s =l .85. Elec- 
tronic energies in CEIMC were calculated using Variational 
Monte Carlo. 



r s 


projection time (a.u.) 1 


time step (a.u.) 


1.05 


0.456 


0.008 


1.10 


0.504 


0.008 


1.15 


0.550 


0.01 


1.25 


0.660 


0.012 


1.40 


0.732 


0.012 


1.55 


0.975 


0.015 



TABLE V: Projection time and time step used in the RQMC 
calculations. 
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